;+
; NAME:
;       HUMAN_VERIFICATION_SLIDES
;
; PURPOSE:
;       Creates the Human Verification Slides needed to corroborate
;       the distance_omnibus algorithmic associations and velocities
;       and everything else.
;
; CATEGORY:
;       BGPS distance-omnibus ancillary routine
;
; CALLING SEQUENCE:
;       human_verification_slides [,csv_file][,/VERBOSE][,POST_DIR=post_dir]
;                                 [,START=start][,REAR=rear][,/ALLONE]
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       CSV_FILE -- Name of the Distance Database CSV file
;                   [Default: 'bgps_distance_database.csv']
;       POST_DIR -- Driectory in which may be found the
;                   postage-stamp images for plotting [Default:
;                   './local/postage/']
;       START    -- Which BGPS source to start with [Default: 1]
;       REAR     -- Which BGPS source to stop with [Default: last]
;       **NOTE**: Maybe also need optional inputs for the various
;                 dependent structures... (presently hardwired)
;
; KEYWORD PARAMETERS:
;       VERBOSE  -- Verbose output
;
; OUTPUTS:
;       Writes .eps slides to the ./local/hv_slides directory, named
;       by Bolocat CNUM (i.e. bgps4313.eps)
;
; OPTIONAL OUTPUTS:
;       NONE
;
; EXAMPLE:
;       human_verification_slides
;
; MODIFICATION HISTORY:
;       Created:  07/09/10, TPEB -- Initial version.
;       Modified: 11/07/10, TPEB -- Make plotted lines thicker.
;       Modified: 01/26/11, TPEB -- Added velocity legend for dense
;                                   gas
;       Modified: 02/01/11, TPEB -- Made conformal changes to use the
;                                   'irdc_full_data.sav' version of
;                                   the Peretto & Fuller (2009) IRDC
;                                   catalog, and added functionality
;                                   to plot multiple IRDCs in panel 5.
;       Modified: 04/27/11, TPEB -- Changed how V_TAN is calculated,
;                                   since it is being compared to
;                                   spectra in V_LSR, not the true
;                                   physical radial velocity (panel 4).
;                                   Also added documentation for
;                                   optional inputs / keywords, and
;                                   some code cleanup.
;       Modified: 07/13/11, TPEB -- Cleaned up the memory hog by
;                                   undefining the big variables at
;                                   the end of the routine.
;       Modified: 09/12/11, TPEB -- Switching to the use of cgColor()
;                                   in the PROBCOLOR.pro calls, and
;                                   conformal changes.
;       Modified: 10/26/11, TPEB -- Added maser location information
;                                   for BGPS sources associated with
;                                   BeSSeL maser parallax measurement.
;                                   Plots a star in panel #1.  Also,
;                                   code cleanup to reflect addition
;                                   of new routines in the
;                                   repository.
;       Modified: 11/08/11, TPEB -- Removing use of the CSO CS(5-4)
;                                   sprectra from use in the
;                                   repository.
;       Modified: 12/09/11, TPEB -- Code cleanup of routines used, and
;                                   moving to cgColor() whenever
;                                   possible.  Removed ALLONE keyword
;
;-

PRO HUMAN_VERIFICATION_SLIDES, csv_file, POST_DIR=postdir, VERBOSE=verbose, $
                               START=start, REAR=rear
  
  ;; Check inputs and keywords
  IF n_elements(csv_file) EQ 0 THEN csv_file = 'bgps_distance_database.csv'
  IF KEYWORD_SET( verbose ) THEN be_verbose = 1b ELSE be_verbose = 0b
  IF n_elements(post_dir) EQ 0 THEN POST_DIR='./local/postage/'
  
  ;; Get repository params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  d = (dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART) / 1.d3
  
  ;; Read in the CSV file
  s = read_bgps_csv(csv_file, csv, VERBOSE=verbose)
  
  ;; Read in the BeSSeL Parallax data
  restore,'./ancillary/bessel_parallax.sav',/ver
  plotsym,3,2,thick=3
  
  ;; **NOTE**: Maybe add optional inputs for these (otherwise
  ;; hardwired) required IDL structures...

  ;; Restore velocity structures
  restore,'./local/bgps_grs_spectra.sav',/ver
  restore,'./local/bgps_hi_spectra.sav',/ver
  restore,'./local/bgps_velocity_struct.sav',/ver
  restore,'./local/bgps_irdc.sav',/ver
  
  ;; Restore probability structures
  restore,'./local/bgps_pvec.sav',/ver
  
  n_prob = N_TAGS(pvec[0]) - 4
  help,pvec[0],/str
  
  ;; Define FOR-loop begining and ending
  IF KEYWORD_SET(start) THEN start=long(start) ELSE start=0L
  IF KEYWORD_SET(rear) THEN rear=long(rear) ELSE rear=csv.nrows-1
  
  FOR i=start, rear DO BEGIN
     
     v = velocity[i]
     titlestr = string(s[i].cnum,s[i].glon_peak,s[i].glat_peak,format=$
                       "('BGPS #',I4,'  l = ',F5.1,'  b = ',F5.2)")
     
     cnum = string(s[i].cnum,format="(I04)")
     
     bgps_post = readfits(POST_DIR+'bgps'+cnum+'.fits',bhdr,/silent)
     bgps_labl = readfits(POST_DIR+'label'+cnum+'.fits',lhdr,/silent)
     bgps_smoo = readfits(POST_DIR+'smooth'+cnum+'.fits',shdr,/silent)
     ;;print,size(bgps_post,/DIM),size(bgps_labl,/DIM),size(bgps_smoo,/DIM)
     
     ;; Check if within GLIMPSE coverage...
     irac_file = POST_DIR+'irac'+cnum+'.fits'
     IF FILE_TEST(irac_file,/READ) THEN BEGIN
        irac_post = readfits(irac_file,ihdr,/silent) 
        have_irac = 1b
     ENDIF ELSE have_irac = 0b 
     igps_file = POST_DIR+'igps'+cnum+'.fits'
     IF FILE_TEST(igps_file,/READ) THEN BEGIN
        igps_post = readfits(igps_file,hihdr,/silent)
        have_igps = 1b
     ENDIF ELSE have_igps = 0b
     
     
     bgps_sz = size(bgps_post, /DIM)
     
     ;; Set Plot Range for BGPS Postage Stamp image
     bmask = bgps_post * bgps_labl
     lmask = bgps_post * (1 - bgps_labl)
     
     IF max(lmask,/nan) GT 2.*max(bmask,/nan) THEN $
        bplot = bmask + lmask / max(lmask,/nan) * 2.*max(bmask,/nan) ELSE $
           bplot = bgps_post
     
     bgps_pr = set_plot_range(bplot)
     
     
     
     ;;=======================================================
     ;; BEGIN PLOTTING
     ;;=======================================================
          
     myps,'local/hv_slides/bgps'+cnum+'.eps',xsize=16.5,ysize=10,ct=3
     print,'Writing local/hv_slides/bgps'+cnum+'.eps'
          
     ;; Set up Multiplot
     erase & multiplot,[3,2],mtitle=titlestr,xgap=0.025,ygap=0.025,$
                       mTitSize=1.5
     
     
     ;;==================================
     ;; Frame 1: BGPS
     plotimage,bgps_post,range=bgps_pr,/preserve_aspect,xst=4,yst=4,min_dpi=60
     
     xc = findgen(bgps_sz[0]) + 0.5
     yc = findgen(bgps_sz[1]) + 0.5
     
     contour,bgps_labl,xc,yc,levels=[0.5],thick=3,$
             color=cgColor('Spring Green'),/overplot
     plot_wcs_axes,bhdr,1,lonarr=blarr,latarr=bbarr
     
     ;; Check against the BeSSeL parallax catalog, and plot a star at
     ;;   the position of associated masers.
     px_ind = WHERE(parallax.cnum EQ s[i].cnum, npx)
     IF npx GT 0 THEN plots,parallax[px_ind].l,parallax[px_ind].b,$
                            psym=8,thick=2,color=cgColor('RYB7')
     
     
     ;;==================================
     ;; Frame 2: GLIMPSE 8-um images (IRAC Band 4)
     multiplot
     IF have_irac THEN BEGIN
        
        ;; Plot GLIMPSE images on a LOG stretch.  (Already checked
        ;;   that there are not enough negative pixels to harm us.)
        plotmap,irac_post,ihdr,AXSEL=1,/log,ct=0
        
        ;;  BGPS flux contours here...
        contour,bgps_labl,blarr,bbarr,levels=[0.5],thick=5,$
                color=cgColor('Spring Green'),/overplot
        cgLoadct,40,/silent
        contour,bgps_smoo,blarr,bbarr,levels=[.025,.05,0.1,0.25,0.5],$
                c_colors=[1,2,3,4,5]*50,/overplot,thick=1.5
        
     ENDIF ELSE BEGIN
        plot,blarr,bbarr,xst=5,yst=5,/nodata,/isotropic
        lpx = 0.14*(!x.crange[1] - !x.crange[0])+!x.crange[0]
        lpy = 0.70*(!y.crange[1] - !y.crange[0])+!y.crange[0]
        al_legend,position=[lpx,lpy],['Outside','GLIMPSE','coverage'],$
                  charthick=2,box=0,charsize=1.5
        plot_wcs_axes,bhdr,1
     ENDELSE
     
     
     ;;==================================
     ;; Frame 3: 21-cm continuum
     multiplot
     IF have_igps THEN BEGIN
        
        ;; Plot IGPS images on a LOG stretch.  (Already checked
        ;; that not enough negative pixels to harm us.)
        suggest_pr = set_plot_range(alog(igps_post),level=0.96)
        igps_pr    = [ALOG(10.),ALOG(!MW.TA_CONT) > suggest_pr[1]]
        plotmap,igps_post,hihdr,AXSEL=1,/log,ct=3,RANGE=igps_pr,$
                LONARR=ilarr,LATARR=ibarr
        cgcolorbar,/vert,/right,charsize=0.65,range=exp(igps_pr),$
                   position=[0.96,0.545,0.97,0.94],format="(F0.1)",$
                   title='21-cm continuum  T!dA!n (K)',/ylog
        
        ;; Add contours at T=100K (continuum limit used by
        ;; Roman-Duaval, et al (2009)
        contour,igps_post,ilarr,ibarr,levels=[ALOG(!MW.TA_CONT)],thick=10,$
                  c_colors=[cgColor('Blue')],/overplot,c_linestyle=[3]
        
        ;;  BGPS source outline and bdrspec region (from HISA.pro) overplots
        contour,bgps_labl,blarr,bbarr,levels=[0.5],thick=5,$
                color=cgColor('Spring Green'),/overplot
        print,size(hi[i].border,/DIM),size(blarr,/DIM),size(bbarr,/DIM)
        contour,hi[i].border,blarr,bbarr,levels=[0.5],thick=2,$
                color=cgColor('Red'),/overplot
        
     ENDIF ELSE BEGIN
        plot,blarr,bbarr,xst=5,yst=5,/nodata,/isotropic
        lpx = 0.14*(!x.crange[1] - !x.crange[0])+!x.crange[0]
        lpy = 0.70*(!y.crange[1] - !y.crange[0])+!y.crange[0]
        al_legend,position=[lpx,lpy],box=0,charsize=1.5,charthick=2,$
                  ['Outside IGPS','21-cm continuum','coverage']
        plot_wcs_axes,bhdr,1
     ENDELSE
     
     
     ;;==================================
     ;; Frame 4: Combined Velocity Spectra
     multiplot
     plot,v.v_std,v.v_std,yr=[-0.25,8.75],/yst,xtit='V!dLSR!n [km s!u-1!n]',$
          ytit='Scaled Intensity',tit=tit,/nodata,ytickname=replicate(' ',30),$
          xr=[-100,200]
     
     black = cgColor('Opposite')
     shift = -1.25
     lcs = 0.75
     vthick = 2
     
     ;; HCOP
     shift += 1.25
     IF TOTAL(v.hcop) NE 0. THEN $
        al_legend,position=[175,shift+0.5],['**'],box=0,charsize=lcs
     scale = max(v.hcop) > 2.d
     oplot,v.v_std,v.hcop  / scale + shift, thick=vthick,color=black
     al_legend,position=[!x.crange[0],shift + 0.5],['HCO!u+!n'],$
               box=0,charsize=lcs
     
     ;; N2HP
     shift += 1.25
     IF TOTAL(v.n2hp) NE 0. THEN $
        al_legend,position=[175,shift+0.5],['**'],box=0,charsize=lcs
     scale = max(v.n2hp) > 2.d
     oplot,v.v_std,v.n2hp  / scale + shift, thick=vthick,color=black
     al_legend,position=[!x.crange[0],shift + 0.5],['N!d2!nH!u+!n'],$
               box=0,charsize=lcs
     
     ;; CS21
     shift += 1.25
     IF TOTAL(v.cs21) NE 0. THEN $
        al_legend,position=[175,shift+0.5],['**'],box=0,charsize=lcs
     scale = max(v.cs21) > 2.d
     oplot,v.v_std,v.cs21  / scale + shift, thick=vthick,color=black
     al_legend,position=[!x.crange[0],shift + 0.5],['CS'],$
               box=0,charsize=lcs
     
     ;; NH3
     shift += 1.25
     IF TOTAL(v.nh3) NE 0. THEN $
        al_legend,position=[175,shift+0.5],['**'],box=0,charsize=lcs
     usemax = 0.5 > max(v.nh3)
     oplot,v.v_std,v.nh3/usemax + shift, thick=vthick,color=black
     al_legend,position=[!x.crange[0],shift + 0.5],['NH!d3!n'],$
               box=0,charsize=lcs
     
     ;; Dense Gas
     shift += 1.25
     yfit = fltarr(n_elements(v.v_std))
     IF TOTAL(v.densegas) NE 0. THEN BEGIN
        max = max(v.densegas)
        yfit = MPFITPEAK(v.v_std, v.densegas, A, NTERMS=3, /POSIT)
        v_densegas = A[1]
        al_legend,position=[175,shift+0.5],['**'],box=0,charsize=lcs
        al_legend,position=[!x.crange[0],shift+0.95],box=0,charsize=lcs,$
                  ['V!dLSR!n = '+string(A[1],format="(F0.1)")+' km s!u-1!n']
     ENDIF ELSE BEGIN
        max=1.
        v_densegas = -1000.
     ENDELSE
     
     ;; oplot,v.v_std,v.densegas / max + shift
     oplot,v.v_std,v.densegas /max + shift, color=probcolor(3), thick=vthick
     al_legend,position=[!x.crange[0],shift + 0.5],['Dense Gas'],$
               box=0,charsize=lcs
     
     ;; GRS 13CO
     shift += 1.25
     IF TOTAL(v.grs13co) NE 0. THEN BEGIN
        max = max(v.grs13co)
        al_legend,position=[175,shift+0.5],['**'],box=0,charsize=lcs
     ENDIF 
     oplot,v.v_std,v.grs13co  / max + shift, color=probcolor(0), thick=vthick
     al_legend,position=[!x.crange[0],shift + 0.5],['GRS !u13!nCO'],$
               box=0,charsize=lcs
     
     ;; HISA
     shift += 1.25
     IF TOTAL(hi[i].spectrum) NE 0 THEN BEGIN
        max = max(hi[i].spectrum)
        al_legend,position=[175,shift+0.5],['**'],box=0,charsize=lcs
     ENDIF ELSE max = 1
     
     oplot,hi[i].v_std,hi[i].spectrum/max + shift,color=probcolor(1), $
           thick=vthick
     al_legend,position=[!x.crange[0],shift + 0.5],['HISA - Invert'],$
               box=0,charsize=lcs
     al_legend,position=[!x.crange[0],shift + 1.0],$
               ['HISA Prob: '+string(velocity[i].p_hisa,format="(F0.3)")],$
               box=0,charsize=lcs
     
     IF v_densegas NE -1000 THEN $
        vline,v_densegas,color=probcolor(3), thick=vthick
     
     ;; Tangent Velocity -- in V_LSR frame, not actual radial velocity
     ;; Calculate VLSR for line-of-sight, and interpolate to tan_dist
     vlsr = vphys2vlsr(s[i].glon_peak, s[i].glat_peak, d*1.d3)
     vtan = INTERPOL( vlsr, d*1.d3, !MW.R0*cos(s[i].glon_peak*!dtor), /LSQUAD)
     
     ;;vtan = !MW.V0*(1. - sin(s[i].glon_peak * !dtor))  ;; OLD
     vline,vtan, thick=vthick
     al_legend,position=[vtan-6,!y.crange[1]*0.98],['V!dtan!n'],$
               box=0,charsize=0.75
     
     
     
     ;;==================================
     ;; Frame 5: IRDC ellipses
     multiplot
     plot,blarr,bbarr,xst=5,yst=5,/nodata,/isotropic
     plot_wcs_axes,bhdr,1
     contour,bgps_labl,blarr,bbarr,levels=[0.5],thick=2,$
             color=cgColor('Spring Green'),/overplot
     
     IF (bgps_irdc[i].n_irdc GT 0) THEN BEGIN
        ;; BGPS Source
        tvellipse,bgps_irdc[i].b_min/3600.,bgps_irdc[i].b_maj/3600.,$
                  s[i].glon_cen,s[i].glat_cen,-bgps_irdc[i].pa_b,$
                  /DATA,color=cgColor('Orange'),thick=5
        
        ;; Loop through how many IRDCs we have for this source
        FOR jj=0L, bgps_irdc[i].n_irdc-1 DO BEGIN
           
           tvellipse,bgps_irdc[i].i_min[jj]/3600.,bgps_irdc[i].i_maj[jj]/3600.,$
                     bgps_irdc[i].i_lcen[jj],bgps_irdc[i].i_bcen[jj],$
                     -bgps_irdc[i].pa_i[jj],$
                     /DATA,color=cgColor('Purple'),thick=5
        ENDFOR
        
        ;; Al_Legend for the ellipse colors
        plsz = [!x.crange[1]-!x.crange[0],!y.crange[1]-!y.crange[0]]
        lp = [0.65*plsz[0] + !x.crange[0],0.15*plsz[1] + !y.crange[0]]
        al_legend,position=lp,['BGPS','IRDC'],color=['Orange','Purple'],$
                  linestyle=[0,0],charsize=0.9,box=0,thick=[5,5]
        
     ENDIF
     plsz = [!x.crange[1]-!x.crange[0],!y.crange[1]-!y.crange[0]]
     lp = [0.05*plsz[0] + !x.crange[0],0.95*plsz[1] + !y.crange[0]]
     al_legend,position=lp,box=0,charsize=1.0,$
               [string(bgps_irdc[i].amplitude,$
                       format="('IRDC Assoc. Prob: ',F0.3)")]
     
     
     
     ;;==================================
     ;; Frame 6: Porbability Functions
     multiplot
     ps = 2b
     plot,d,d*0.+.025,/NODATA,$
          xtitle='Distance (kpc)',ytitle='Probability per 20 pc bin',$
          ytickformat="(F0.3)",color=cgColor('Opposite')
     al_legend,position=[0.60*!x.crange[1],0.95*!y.crange[1]],$
               ['GRS !u13!nCO','HISA','IRDC','KDIST','H!d2!n','MASER','Total'],$
               color=probcolor( [0,1,2,3,4,5,-100], /NAME ), $
               linestyle=0,thick=[4,4,4,4,4,4,5], number=0.75
     al_legend,position=[0.60*!x.crange[1],0.60*!y.crange[1]],$
               [string(pvec[i].disc.nearp,format="('P!dnear!n = ',F5.3)")],$
               box=0
     
     FOR kk=1L, n_prob DO $
        oplot,d,pvec[i].(kk),color=probcolor(kk-1),thick=4
     
     oplot,d,pvec[i].total,color=probcolor(-100),thick=5
     
     tandist = !MW.R0*cos(s[i].glon_peak * !dtor)/1.d3
     ;; Only plot if within box
     IF (tandist GE !x.crange[0]) AND (tandist LE !x.crange[1]) THEN BEGIN
        plots,tandist,!y.crange[0]
        plots,tandist,!y.crange[1],/continue,linestyle=2
     ENDIF
     
     
     ;; Reset for next object
     multiplot,/reset
     
     myps,/done
     
  ENDFOR
  
  ;; Clean up the memory, dude
  undefine,pvec
  undefine,velocity
  undefine,bgps_irdc
  undefine,grs
  undefine,hi
  
END
